home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
answrbok
/
6_10.lha
/
6_10
/
6_10div.c
< prev
next >
Wrap
Text File
|
1993-08-08
|
14KB
|
534 lines
* Copyright (c) 1990 by AT&T Bell Telephone Laboratories, Incorporated. */
* The C++ Answer Book */
* Tony Hansen */
* All rights reserved. */
*
Divide u[1..m+n] by v[1..n] to
form q[0..m] and r[1..n]
Based on:
The Art of Computer Programming, volume 2
D. Knuth, Section 4.3.1, Algorithm D and
exercise 16.
/
include <lint.h>
include <debug.h> /* DELETE */
oid dodivmod(LINT u, LINT v,
LINT "ient, LINT &remainder)
int negquotient = 0, negremainder = 0;
{ if (debug > 1) cerr << "\n"; /* DELETE */ }
{ if (debug > 1) cerr << form("u=%4.4x%4.4x%4.4x%4.4x", u.s[0], u.s[1], u.s[2], u.s[3]) << "\n"; /* DELETE */ }
{ if (debug > 1) cerr << form("v=%4.4x%4.4x%4.4x%4.4x", v.s[0], v.s[1], v.s[2], v.s[3]) << "\n"; /* DELETE */ }
/*
Make u (the dividend) positive.
The sign of the remainder will
match the sign of the dividend.
*/
if (u.isneg())
{
negremainder = 1;
u = -u;
if (debug) cerr << "PT 1\n"; /* DELETE */ }
{ if (debug > 1) cerr << "neg remainder\n"; /* DELETE */ }
}
else /* DELETE */
{ /* DELETE */
if (debug) cerr << "PT 2\n"; /* DELETE */ }
} /* DELETE */
/*
Make v (the divisor) positive.
The sign of the quotient will be
negative if the sign of the divisor
and dividend do not match, else
positive.
*/
if (v.isneg())
{
if (debug) cerr << "PT 3\n"; /* DELETE */ }
negquotient = !negremainder;
v = -v;
}
else
{ /* DELETE */
if (debug) cerr << "PT 4\n"; /* DELETE */ }
negquotient = negremainder;
} /* DELETE */
{ if (debug > 1) cerr << (negquotient ? "negquotient\n" : ""); /* DELETE */ }
// set local variables
for (int m_n = 4, uoffset = 0;
uoffset < 4 && u.s[uoffset] == 0;
m_n--, uoffset++)
{ /* DELETE */
if (debug) cerr << "PT 5\n"; /* DELETE */ }
;
} /* DELETE */
for (int n = 4, voffset = 0;
voffset < 4 && v.s[voffset] == 0;
n--, voffset++)
{ /* DELETE */
if (debug) cerr << "PT 6\n"; /* DELETE */ }
;
} /* DELETE */
int m = m_n - n;
{ if (debug > 1) cerr << "m=" << m <<", n=" << n << ", m_n=" << m_n << "\n"; /* DELETE */ }
/*
If n == 0, then division by zero
*/
if (n == 0)
{
if (debug) cerr << "PT 7\n"; /* DELETE */ }
{ if (debug > 1) cerr << "division by zero\n"; /* DELETE */ }
error("division by zero!");
quotient = remainder = u;
return;
}
/*
For n == 1, use simpler algorithm
*/
else if (n == 1)
{
if (debug) cerr << "PT 8\n"; /* DELETE */ }
{ if (debug > 1) cerr << "n==1\n"; /* DELETE */ }
// See Exercise 16 after Section 4.3.1
LINT q = 0;
LINT_Ltype prevu = 0;
LINT_type v1 = v.s[3];
for (int r = uoffset; r < m_n + uoffset; r++)
{
if (debug) cerr << "PT 9\n"; /* DELETE */ }
LINT_Ltype t = u.s[r] + prevu * LINT_base;
LINT_type tmpq = LINT_type (t / v1);
q.s[r] = tmpq; // % LINT_base
prevu = t - v1 * tmpq;
}
if (negquotient)
{ /* DELETE */
if (debug) cerr << "PT 10\n"; /* DELETE */ }
quotient = -q;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 11\n"; /* DELETE */ }
quotient = q;
} /* DELETE */
if (negremainder)
{ /* DELETE */
if (debug) cerr << "PT 12\n"; /* DELETE */ }
remainder = -prevu;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 13\n"; /* DELETE */ }
remainder = prevu;
} /* DELETE */
return;
}
/*
Degenerate case of length(u) < length(v)
i.e., m < 0, implying that u < v
*/
else if (m_n < n)
{
if (debug) cerr << "PT 14\n"; /* DELETE */ }
{ if (debug > 1) cerr << "m_n < n\n"; /* DELETE */ }
quotient = 0L;
if (negremainder)
{ /* DELETE */
if (debug) cerr << "PT 15\n"; /* DELETE */ }
remainder = -u;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 16\n"; /* DELETE */ }
remainder = u;
} /* DELETE */
return;
}
/*
Degenerate case of length(u) == length(v)
i.e., m == 0, possibly implying that
u < v or u == v
*/
else if (m_n == n)
{
if (debug) cerr << "PT 17\n"; /* DELETE */ }
{ if (debug > 1) cerr << "m_n == n\n"; /* DELETE */ }
int cmp = LINT_cmp(&u.s[uoffset], &v.s[voffset], m_n);
if (cmp < 0) // u < v
{
if (debug) cerr << "PT 18\n"; /* DELETE */ }
{ if (debug > 1) cerr << "cmp < 0\n"; /* DELETE */ }
quotient = 0L;
if (negremainder)
{ /* DELETE */
if (debug) cerr << "PT 19\n"; /* DELETE */ }
remainder = -u;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 20\n"; /* DELETE */ }
remainder = u;
} /* DELETE */
return;
}
else if (cmp == 0) // u == v
{
if (debug) cerr << "PT 21\n"; /* DELETE */ }
{ if (debug > 1) cerr << "cmp == 0\n"; /* DELETE */ }
if (negquotient)
{ /* DELETE */
if (debug) cerr << "PT 22\n"; /* DELETE */ }
quotient = -1L;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 23\n"; /* DELETE */ }
quotient = 1L;
} /* DELETE */
remainder = 0L;
return;
}
else /* DELETE */
{ /* DELETE */
if (debug) cerr << "PT 24\n"; /* DELETE */ }
} /* DELETE */
}
else /* DELETE */
{ /* DELETE */
if (debug) cerr << "PT 25\n"; /* DELETE */ }
} /* DELETE */
/*
Now call out all of the guns from Knuth
*/
/*
D1(a) [Normalize.]
Set d <- b/(v1+1).
*/
LINT_type d = LINT_base / (v.s[voffset] + 1);
{ if (debug > 1) cerr << "d=" << form("0x%4.4x", d) << "\n"; /* DELETE */ }
/*
D1(b) [Normalize.]
Set (u[0]u[1]...u[m+n]) to
(u[1]u[2]...u[m+n])*d
*/
LINT_type uv[5];
LINT_type k;
if (d == 1)
{
if (debug) cerr << "PT 26\n"; /* DELETE */ }
{ if (debug > 1) cerr << "d == 1, copy old u value\n"; /* DELETE */ }
// copy old value
uv[0] = 0;
memcpy((char*)&uv[1], (char*)&u.s[uoffset],
m_n * sizeof(LINT_type));
}
else
{
if (debug) cerr << "PT 27\n"; /* DELETE */ }
{ if (debug > 1) cerr << "multiply u by d\n"; /* DELETE */ }
// multiply u by d
k = 0;
for (int Oi = m_n - 1 + uoffset, i = m_n;
i > 0; Oi--, i--)
{
if (debug) cerr << "PT 28\n"; /* DELETE */ }
LINT_Ltype t = u.s[Oi] * d + k;
uv[i] = t; // % LINT_base
k = t / LINT_base;
{ if (debug > 1) cerr << "done setting u.s[" << Oi << "], uv[" << i << "] = " << form("0x%4.4x", uv[i]) << "\n"; /* DELETE */ }
}
uv[i] = k;
}
{ if (debug > 1) { cerr << "uv=0x"; for (int jj=0; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
/*
D1(c) [Normalize.]
Set (v[1]v[2]...v[n]) to
(v[1]v[2]...v[n]) * d
*/
LINT_type vv[4];
if (d == 1)
{
if (debug) cerr << "PT 29\n"; /* DELETE */ }
{ if (debug > 1) cerr << "d == 1, copy old v value\n"; /* DELETE */ }
// copy old value
memcpy((char*)&vv[0], (char*)&v.s[voffset],
n * sizeof(LINT_type));
}
else
{
if (debug) cerr << "PT 30\n"; /* DELETE */ }
{ if (debug > 1) cerr << "multiply v by d\n"; /* DELETE */ }
// multiply v by d
k = 0;
for (int Oi = n - 1 + voffset, i = n - 1;
i >= 0; Oi--, i--)
{
if (debug) cerr << "PT 31\n"; /* DELETE */ }
LINT_Ltype t = v.s[Oi] * d + k;
vv[i] = t; // % LINT_base;
k = t / LINT_base;
{ if (debug > 1) cerr << "done setting v.s[" << Oi << "], vv[" << i << "] = " << form("0x%4.4x", vv[i]) << "\n"; /* DELETE */ }
}
}
{ if (debug > 1) { cerr << "vv=0x"; for (int jj=0; jj < n; jj++) cerr << form("%4.4x", vv[jj]); cerr << "\n"; } /* DELETE */ }
/*
D2 [Initialize j]
Set j <- 0
D7 [Loop on j]
Set j <- j+1
Loop if j <= m
*/
LINT_type v1 = vv[0];
LINT_type v2 = vv[1];
{ if (debug > 1) cerr << "v1=" << form("0x%4.4x", v1) << ", v2=" << form("0x%4.4x", v2) << "\n"; /* DELETE */ }
LINT_type qv[5];
for (int j = 0; j <= m; j++)
{
if (debug) cerr << "PT 32\n"; /* DELETE */ }
{ if (debug > 1) cerr << "j=" << j << "\n"; /* DELETE */ }
/*
D3 [Calculate q^]
If u[j] == v[1]
Set q^ <- base - 1
Else
Set q^ <-
(u[j] * base + u[j+1]) / v[1]
While v[2] * q^ >
(u[j] * base + u[j+1] - q^ * v[1]) *
base + u[j+2]
Set q^ <- q^ - 1
*/
LINT_Ltype q_hat;
if (uv[j] == v1)
{ /* DELETE */
if (debug) cerr << "PT 33\n"; /* DELETE */ }
q_hat = LINT_base - 1;
{ if (debug > 1) cerr << "uv[" << j << "]=" << form("0x%4.4x", uv[j]) << " == v1, q^=" << form("0x%4.4x", q_hat) << "\n"; /* DELETE */ }
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 34\n"; /* DELETE */ }
q_hat = (uv[j] * LINT_base + uv[j+1]) / v1;
{ if (debug > 1) cerr << "uv[" << j << "]=" << form("0x%4.4x", uv[j]) << " != v1, q^=" << form("0x%4.4x", q_hat) << "\n"; /* DELETE */ }
} /* DELETE */
LINT_Ltype u_j = uv[j];
LINT_Ltype u_j1 = uv[j+1];
LINT_Ltype u_j2 = uv[j+2];
for ( ; ; q_hat--)
{
if (debug > 1) cerr << "q^=" << form("0x%4.4x", q_hat) << "\n"; /* DELETE */ }
/*
The following statements build up the
following test which combines with the
for loop above to give the while
statement mentioned in the last comment.
if ((v2 * q_hat) <=
((u_j * LINT_base + u_j1 - v1 *
q_hat) * LINT_base + u_j2))
q^--;
*/
// u_j_q_hat = u_j * LINT_base + u_j1
union {
LINT_Ltype L[2];
LINT_type S[4];
} u_j_q_hat;
u_j_q_hat.L[0] = 0;
u_j_q_hat.S[2] = u_j;
u_j_q_hat.S[3] = u_j1;
// ... - v1 * q_hat
u_j_q_hat.L[1] -= v1 * q_hat;
if (u_j_q_hat.S[2] != 0)
break;
// ... * LINT_base
u_j_q_hat.S[2] = u_j_q_hat.S[3];
// ... + u_j2
u_j_q_hat.S[2] = u_j2;
if ((v2 * q_hat) <= u_j_q_hat.L[1])
break;
if (debug) cerr << "PT 35\n"; /* DELETE */ }
if (debug > 1) cerr << "q^--\n"; /* DELETE */ }
}
/*
D4 [Multiply and subtract.]
Replace u[j..j+n] by
u[j..j+n] -
q^ * v[1..n]
*/
// set nv <- q^ * (v[1..n])
LINT_type nv[5];
k = 0;
for (int dl = n, vl = n-1; vl >= 0; vl--, dl--)
{
if (debug) cerr << "PT 36\n"; /* DELETE */ }
LINT_Ltype t = vv[vl] * q_hat + k;
nv[dl] = t; // % LINT_base;
k = t / LINT_base;
{ if (debug > 1) cerr << "done setting nv[" << dl << "] = " << form("0x%4.4x", nv[dl]) << ", vv[" << vl << "] = " << form("0x%4.4x", vv[vl]) << "\n"; /* DELETE */ }
}
nv[0] = k;
{ if (debug > 1) cerr << "done setting nv[" << dl << "] = " << form("0x%4.4x", nv[dl]) << "\n"; /* DELETE */ }
{ if (debug > 1) { cerr << "nv=0x"; for (int jj=0; jj <= n; jj++) cerr << form("%4.4x", nv[jj]); cerr << "\n"; } /* DELETE */ }
// subtract nv[0..n] from u[j..j+n]
int borrow = 0;
int ul = j + n;
{ if (debug > 1) cerr << "ul=" << ul << ", dl=" << n << "\n"; /* DELETE */ }
for (dl = n; dl >= 0; dl--, ul--)
{
if (debug) cerr << "PT 37\n"; /* DELETE */ }
LINT_Ltype t = uv[ul] - nv[dl] - borrow;
uv[ul] = t; // % LINT_base
borrow = (t / LINT_base) ? 1 : 0;
{ if (debug > 1) cerr << "done setting uv[" << ul << "] = " << form("0x%4.4x", uv[ul]) << "\n"; /* DELETE */ }
}
{ if (debug > 1) { cerr << "uv=0x"; for (int jj=0; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
{ if (debug > 1) cerr << "borrow=" << borrow << "\n"; /* DELETE */ }
/*
D5 [Test remainder]
Set q[j] <- q^
if u[j] < 0
D6 [Add back]
add 0v[1..n] back to u[j..j+n]
q[j]--
*/
qv[j] = q_hat;
{ if (debug > 1) cerr << "qv[" << j << "] = " << form("0x%4.4x", qv[j]) << "\n"; /* DELETE */ }
if (borrow != 0)
{
if (debug) cerr << "PT 38\n"; /* DELETE */ }
{ if (debug > 1) cerr << "borrow != 0, add v back to u\n"; /* DELETE */ }
for (k = 0, ul = j + n, vl = n - 1;
vl >= 0;
vl--, ul--)
{
if (debug) cerr << "PT 39\n"; /* DELETE */ }
LINT_Ltype t = uv[ul] + vv[vl] + k;
uv[ul] = t; // % LINT_base
{ if (debug > 1) cerr << "done setting uv[" << ul << "] = " << form("0x%4.4x", uv[ul]) << "\n"; /* DELETE */ }
k = t / LINT_base;
}
uv[j] += k;
{ if (debug > 1) cerr << "done setting uv[" << ul << "] = " << form("0x%4.4x", uv[ul]) << "\n"; /* DELETE */ }
{ if (debug > 1) { cerr << "uv=0x"; for (int jj=0; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
qv[j]--;
{ if (debug > 1) cerr << "qv[" << j << "] = " << form("0x%4.4x", qv[j]) << "\n"; /* DELETE */ }
}
else /* DELETE */
{ /* DELETE */
if (debug) cerr << "PT 40\n"; /* DELETE */ }
} /* DELETE */
}
/*
D8 [Unnormalize]
q[0..m] is quotient
u[m+1..m+n] / d is remainder
*/
LINT q = 0;
memcpy((char*)&q.s[3-m], (char*)&qv[0],
(m+1) * sizeof(LINT_type));
{ if (debug > 1) cerr << "q = " << q << "\n"; /* DELETE */ }
if (negquotient)
{ /* DELETE */
if (debug) cerr << "PT 41\n"; /* DELETE */ }
quotient = -q;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 42\n"; /* DELETE */ }
quotient = q;
} /* DELETE */
// divide u[m+1..m+n] by d
LINT r = 0;
{ if (debug > 1) { cerr << "uv[m+1..m+n]=0x"; for (int jj=m+1; jj <= m_n; jj++) cerr << form("%4.4x", uv[jj]); cerr << "\n"; } /* DELETE */ }
if (d == 1) // nothing special to do
{ /* DELETE */
if (debug) cerr << "PT 43\n"; /* DELETE */ }
memcpy((char*)&r.s[4 - n], (char*)&uv[m + 1],
n * sizeof(LINT_type));
} /* DELETE */
else
{
if (debug) cerr << "PT 44\n"; /* DELETE */ }
LINT_Ltype prevu = 0;
// do division by single digit
for (int rl = 4 - n, ul = m + 1;
ul <= m_n;
ul++, rl++)
{
if (debug) cerr << "PT 45\n"; /* DELETE */ }
LINT_Ltype t = uv[ul] + prevu * LINT_base;
LINT_type tmpq = t / d;
r.s[rl] = tmpq; // % LINT_base
prevu = t - d * tmpq;
}
}
{ if (debug > 1) cerr << "r = " << r << "\n"; /* DELETE */ }
if (negremainder)
{ /* DELETE */
if (debug) cerr << "PT 46\n"; /* DELETE */ }
remainder = -r;
} /* DELETE */
else
{ /* DELETE */
if (debug) cerr << "PT 47\n"; /* DELETE */ }
remainder = r;
} /* DELETE */
INT operator%(LINT u, LINT v)
LINT quot, rem;
dodivmod(u, v, quot, rem);
return rem;
INT operator/(LINT u, LINT v)
LINT quot, rem;
dodivmod(u, v, quot, rem);
return quot;